Lab Session 1

retaildata <- read.csv("retail.csv")
mytimeseries <- ts(retaildata[,4], frequency=12, start=c(1982,4))
autoplot(mytimeseries)

Lab Session 2

ggseasonplot(mytimeseries)

ggsubseriesplot(mytimeseries)

There is a strong trend to around 2008, weak seasonality (with a jump in December each year), and some evidence of cycles (with dips in 2001 and 2012).

Lab Session 3

autoplot(bicoal)

gglagplot(bicoal)

ggAcf(bicoal)

autoplot(dole)

ggseasonplot(dole)

ggsubseriesplot(dole)

gglagplot(dole)

ggAcf(dole)

Lab Session 4

beer <- window(ausbeer, start=1992)
fc <- snaive(beer)
autoplot(fc)

res <- residuals(fc)
autoplot(res)

checkresiduals(fc)

## 
##  Ljung-Box test
## 
## data:  Residuals from Seasonal naive method
## Q* = 32.269, df = 8, p-value = 8.336e-05
## 
## Model df: 0.   Total lags used: 8

There is some remaining autocorrelation in the residuals: the Null of no joint autocorrelation is clearly rejected. We can also see a significant spike on the seasonal (4th lag) in the ACF. There is considerable information remaining in the residuals which has not been captured with the seasonal naïve method. The residuals do not appear to be too far from Normally distributed.

Lab Session 5

train <- window(mytimeseries, end=c(2010,12))
test <- window(mytimeseries, start=2011)
autoplot(cbind(Training=train,Test=test))

f1 <- snaive(train, h=length(test))
autoplot(f1) + autolayer(test)

accuracy(f1, test)

##                    ME     RMSE      MAE      MPE      MAPE      MASE
## Training set 5.502402 27.74935 19.33784 3.369450 10.447161 1.0000000
## Test set     5.227586 23.40458 18.60000 1.619239  8.172045 0.9618449
##                   ACF1 Theil's U
## Training set 0.8703252        NA
## Test set     0.4544630 0.9518744

The number to look at here is the test set RMSE. That provides a benchmark for comparison when we try other models.

checkresiduals(f1)

## 
##  Ljung-Box test
## 
## data:  Residuals from Seasonal naive method
## Q* = 1256.8, df = 24, p-value < 2.2e-16
## 
## Model df: 0.   Total lags used: 24

The residuals do not look like white noise there are lots of dynamics left over that need to be explored. They also do not look close to normal, with very long tails.

The accuracy measure are always sensitive to the training/test split. There are better ways to check the robustness of the methods in terms of accuracy such as using a tsCV().

Lab Session 6

e_mean <- tsCV(mytimeseries, meanf)
e_naive <- tsCV(mytimeseries, naive)
e_snaive <- tsCV(mytimeseries, snaive)
e_drift <- tsCV(mytimeseries, rwf, drift=TRUE)

# Construct squared CV errors matrix
e2 <- cbind(Mean=e_mean^2, Naive=e_naive^2, 
        SNaive=e_snaive^2, Drift=e_drift^2)
# Remove rows with any missing for a fair comparison
e2 <- na.omit(e2)
# Find MSE values
colMeans(e2)

##      Mean     Naive    SNaive     Drift 
## 4685.8665  410.0438  691.7144  411.4290

Lab Session 7

# Best model:
ses(eggs, h=100) %>% accuracy()

##                     ME     RMSE      MAE       MPE     MAPE      MASE
## Training set -2.740622 26.56395 19.38206 -2.872958 10.05648 0.9560879
##                      ACF1
## Training set -0.005983602

holt(eggs, h=100, damped=FALSE) %>% accuracy()

##                      ME     RMSE      MAE       MPE     MAPE      MASE
## Training set 0.04499087 26.58219 19.18491 -1.142201 9.653791 0.9463626
##                    ACF1
## Training set 0.01348202

holt(eggs, h=100, damped=TRUE) %>% accuracy()

##                     ME     RMSE     MAE       MPE     MAPE      MASE
## Training set -2.891496 26.54019 19.2795 -2.907633 10.01894 0.9510287
##                      ACF1
## Training set -0.003195358

These RMSE values are not really comparable because the models have different numbers of parameters.

The best model is the last one.

holt(eggs, h=100, damped=TRUE) %>% forecast() %>% autoplot()

holt(eggs, h=100, damped=TRUE) %>% checkresiduals()

## 
##  Ljung-Box test
## 
## data:  Residuals from Damped Holt's method
## Q* = 10.782, df = 5, p-value = 0.05587
## 
## Model df: 5.   Total lags used: 10

The residuals are pretty good and (just) pass the Ljung-Box test. However, the forecast make little sense (with the prediction intervals going negative very quickly).

Lab Session 8

autoplot(mytimeseries)

The seasonal variation increases with the level of the series. So we need to use multiplicative seasonality.

fc1 <- hw(mytimeseries, seasonal='multiplicative', damped=FALSE)
autoplot(fc1)

fc2 <- hw(mytimeseries, seasonal='multiplicative', damped=TRUE)
autoplot(fc2)

accuracy(fc1)

##                      ME     RMSE      MAE        MPE     MAPE      MASE
## Training set -0.2054187 9.280748 6.244629 -0.2790438 3.338901 0.3336358
##                    ACF1
## Training set 0.05182697

accuracy(fc2)

##                     ME     RMSE      MAE       MPE     MAPE      MASE
## Training set 0.5050393 9.194427 6.128123 0.2458937 3.290157 0.3274111
##                     ACF1
## Training set -0.04640084

There is not much difference between these models, and we would expect the trend to continue, so I would prefer to use the non-damped version.

checkresiduals(fc1)

## 
##  Ljung-Box test
## 
## data:  Residuals from Holt-Winters' multiplicative method
## Q* = 41.483, df = 8, p-value = 1.693e-06
## 
## Model df: 16.   Total lags used: 24

There are significant correlations in the residuals, especially at lag 12. So these residuals do not look like white noise.

train %>% window(end=c(2010,12)) %>%
  hw(seasonal='multiplicative', damped=FALSE) %>%
  accuracy(x=mytimeseries)

##                      ME      RMSE       MAE         MPE     MAPE      MASE
## Training set -0.1212843  9.260807  6.054543 -0.03697659 3.442954 0.3130931
## Test set     10.1364618 19.695875 15.215508  4.30118709 6.977221 0.7868257
##                    ACF1 Theil's U
## Training set 0.09051818        NA
## Test set     0.61161198  0.909366

The test set RMSE is much better than for the seasonal naïve method.

Lab Session 9

fit <- ets(mytimeseries)
summary(fit)

## ETS(M,Ad,M) 
## 
## Call:
##  ets(y = mytimeseries) 
## 
##   Smoothing parameters:
##     alpha = 0.7718 
##     beta  = 0.0025 
##     gamma = 0.0561 
##     phi   = 0.98 
## 
##   Initial states:
##     l = 64.0899 
##     b = 0.2813 
##     s = 0.9869 0.9362 1.0043 1.1705 1.0216 1.0293
##            0.9848 0.9998 0.993 0.9436 0.9694 0.9606
## 
##   sigma:  0.0464
## 
##      AIC     AICc      BIC 
## 4404.041 4405.697 4477.272 
## 
## Training set error measures:
##                     ME     RMSE      MAE       MPE     MAPE     MASE
## Training set 0.5339249 9.613425 6.241919 0.1970991 3.257992 0.333491
##                     ACF1
## Training set -0.03180379

This is equivalent to a damped-trend multiplicative Holt-Winters’ method. The very small \(\beta\) and \(\gamma\) values show that the seasonality and trend change slowly. There is also very little damping (\(\phi\) is close to 1).

bicoal %>% ets() %>% forecast() %>% autoplot()

Not too bad given the noisy data.

bricksq %>% ets() %>% forecast() %>% autoplot()

These are great forecasts which have adapted to the trend and seasonality really well.

dole %>% ets() %>% forecast() %>% autoplot()

The seasonality seems to be over-done, but the wide prediction intervals are probably sensible apart from the fact that they go negative.

lynx %>% ets() %>% forecast() %>% autoplot()

These data are cyclic, and ETS does not handle cycles in time series.